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The random-phase approximation to the ground state correlation energy (RPA) in combination 
with exact exchange (EX) has brought Kohn-Sham (KS) density functional theory one step closer 
towards a universal, "general purpose first principles method" . In an effort to systematically assess 
the influence of several correlation energy contributions beyond RPA, this work presents dissociation 
energies of small molecules and solids, activation energies for hydrogen transfer and non-hydrogen 
transfer reactions, as well as reaction energies for a number of common test sets. We benchmark 
EX+RPA and several flavors of energy functionals going beyond it: second-order screened exchange 
(SOSEX), single excitation (SE) corrections, renormalized single excitation (rSE) corrections, as 
well as their combinations. Both the single excitation correction as well as the SOSEX contribution 
to the correlation energy significantly improve upon the notorious tendency of EX -f RPA to under- 
bind. Surprisingly, activation energies obtained using EX-I-RPA based on a KS reference alone are 
remarkably accurate. RPA-|-SOSEX-|-rSE provides an equal level of accuracy for reaction as well as 
activation energies and overall gives the most balanced performance, which makes it applicable to 
a wide range of systems and chemical reactions. 



I. INTRODUCTION 

In the context of first principles electronic structure 
theory "exact-exchange plus correlation in the random- 
phase approximation (EX+RPA)"— has recently gen- 
erated renewed and widespread interest.—^— In practice, 
the RPA calculations are most often carried out in a non- 
self-consistent manner where the exchange-correlation 
(xc) energy contributions are evaluated with input or- 
bitals corresponding to an approximate, usually semilo- 
cal xc energy functional. The great interest in EX -f RPA 
is largely due to its three attractive features: (z) The 
exact-exchange energy (EX) cancels the spurious self- 
interaction error present in the Hartree energy, (ii) the 
RPA correlation energy is fully non-local and includes 
long-range van dor Waals (vdW) interactions automat- 
ically and seamlessly^ and (Hi) EX-I-RPA is applica- 
ble to small-gap or metallic systems by summing up the 
sequence of "ring" diagrams to infinite order. The lat- 
ter is in contrast to order-by-order perturbation theo- 
ries [e.g. 2nd-order M0ller-Plesset (MP2)^] which break 
down for systems with zero gap. Moreover one can in- 
terpret the RPA as an approach that screens the non- 
local exchange resulting in a frequency dependent non- 
local screened exchange interaction, as opposed to con- 
ventional or global hybrid functionals where the param- 
eters that reduce or "screen" the exact-exchange con- 
tribution are fixed and system independent<^"— Such a 
system independent "screening" is expected to be unre- 
liable for metals or wide gap insulators, where non-local 
exchange is almost entirely screened (metals) or prevails 



to a large extent (insulators). 

While a critical assessment of EX-I-RPA is 
emerging^"— some shortcomings have been 
known for a while. Total energies are typically signif- 
icantly overestimated r^'iii^'^i^ which is caused by 
an overestimation of the correlation energy at the short 
range. Binding energies, on the other hand, show a 
tendency to be underestimated4iii^'i^ii^'i^iiii2^i^ More- 
over, the RPA correlation energy is not self-correlation 
free.^^^ 

It has been demonstrated that the overestimation 
of the absolute correlation energy can be almost en- 
tirely removed by adding a second-order screened ex- 
change (SOSEX) termi^i^i^ For one-electron systems 
the self-interaction error in EX-I-RPA is exactly can- 
celed by adding this tcrm^^i^ however, for systems with 
more than one electron, a many-electron self-interaction 
prevails.— The SOSEX can also be interpreted 
as a correction to the RPA correlation energy that can 
be included to approximately restore the antisymmetry 
of the many electron description)^ Furthermore, SOSEX 
improves binding energies, although a sizeable underesti- 
mation persistsi^i^'^i^'^i^ The underbinding problem 
can also be alleviated, in particular for weakly interact- 
ing systems, by adding a correction deriving from sin- 
gle excitations (SE)^ to EX+RPA built on a reference 
state obtained from Kohn-Sham (KS) density-functional 
theory (DFT). This suggests that RPA(+SOSEX) yields 
good estimations for the correlation energy but errors in 
the exchange energy are sizeable if Kohn-Sham orbitals 
are used to evaluate the exact exchange. 
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In light of these observations it is timely to extend 
the critical assessment of EX+RPA to a wider class of 
systems and to consider combinations of the corrections 
suggested before. In this paper we will address this ob- 
jective by performing benchmark calculations for atom- 
ization energies on an appreciable test set of archety- 
pal insulating solids and small moleculesr^"— as well as 
reaction and activation energies for hydrogen and non- 
hydrogen transfer reactions 4^1^ The schemes we include 
are EX-I-RPA based on KS-DFT reference states, and 
those beyond EX-I-RPA by adding corrections from SE 
or SOSEX individually, or both of them. In addition we 
will also assess the hybrid-type schemes^^ where one re- 
places the total energy at the EX level evaluated with 
KS-DFT orbitals by that evaluated with Hartree-Fock 
(HE) orbitals, as an effective way to approximate the SE 
contribution.— The second-order single excitation cor- 
rection can diverge when the gap between occupied and 
virtual states closes, with detrimental effects for the de- 
scription of the transition states in chemical reactions. As 
briefly discussed in Ref. [2^, including higher-order terms 
in the spirit of RPA permits a resummation of the SE 
correction, as will be demonstrated in Section Hi CI This 
so called renormalized SE (rSE) is well behaved and is 
included in our benchmark tests. 

The paper is organized as follows: Sections |ll] and IIIII 
briefly summarize the important aspects of the under- 
lying theory and the computational parameters of our 
work. Results on molecular and solid-state atomization 
energies as well as reaction energies and barrier heights 
are presented in Section IIVI before we draw conclusions 
m Section El 



II. THEORY 

A. Basics on RPA 

In order to properly position the methods applied in 
the present work within the formal framework of DFT, 
we briefly recapitulate essential equations and outline the 
structure of the functionals used. Currently, for total 
energy calculations, RPA-based functionals usually use 
either KS-DFT reference states, i.e. single-particle wave- 
functions and eigenvalues or generalized KS (GKS)^ ref- 
erence states to compute the nonlocal EX energy as well 
as the nonlocal correlation energyj^i^ii Within this con- 
text, the total energy is defined as 

E[n] = T,[{^,}]+E^[n]+E,^t[n]+E^[{(^,}]+E,[{(j),}], 

(1) 

where the terms deriving from the potential contribu- 
tions in the Hamiltonian, i^n, the electrostatic Hartree 
or Coulomb energy and i?cxt, the (external) electron-ion 
interaction depend on the local density, whereas the last 
two terms, EX energy Ey^ and correlation energy E^, are 
nonlocal contributions. Note that the KS kinetic en- 
ergy, in analogy to the exact-exchange energy, is not 



an explicit functional of the density, but rather of the 
KS orbitals. The nonlocality in E^ is due to the nonlo- 
cal exchange operator acting on each (occupied) orbital 
ipiai'r) associated to spin a and its well known depen- 
dence on the (nonlocal) reduced one-particle density ma- 
trix p^(r,r') = YIT" ^^''^'^'^ 



Ey 



Ip<t(i", I"') 



(2) 



In contrast to the optimized effective potential (OEP) 
method;^"— in HE theory the exchange operator is fully 
nonlocal, and the action of the exchange operator on a 
single-particle wavefunction (i.e. orbital) depends on the 
value of that very orbital throughout the entire space (see 
Ref. Issh . Note that the correlation energy E^ is a func- 
tional of both occupied as well as unoccupied eigenstates 
and requires knowledge of the associated eigenenergies as 
well (see below). However, both, i?x and Ec are implicit 
functionals of the electron density n (see e.g. Ref. [H^. 
Recent work pursuing the construction of a local RPA 
correlation potential are presented in Refs. 57—62. Work 
in this direction is of great value, since it ultimately en- 
ables calculations of a self-consistent RPA correlation en- 
ergies staying rigorously within the KS-DFT picture. 

The RPA correlation energy can be conveniently 
derived from (i) perturbation theory or (ii) from 
the adiabatic-connection fluctuation-dissipation (ACED) 
theorem)^"— Fundamental to the formalism is the adi- 
abatic connection between the Hamiltonian H of an in- 
teracting many-electron system and the corresponding 
noninteracting KS Hamiltonian Hks- Formally, both 
systems may be simultaneously described by a coupling 
constant dependent Hamiltonian H{X) with A being the 
coupling-constant or the scaling factor in the electron- 
electron interaction, vx = Xv{r — r'). The electrons move 
in a A-dependent external potential v^^^. (r). Note that 
the ground-state density of H{X) for all A £ [0, 1] is 
constant and equals the physical ground-state density 
n{r), i.e. the ground-state density of the real system. 
H{X = 1) is the physical many-electron Hamiltonian with 
v\=i{r) ~ Uext(r), and H{X = 0) is the KS Hamiltonian 
with vx=o{r) = vks(i") = Woxt(r) + WH(r) + Vxc(r)- WH(r) 
is the electrostatic Hartree potential and Uxc(r) is the xc 
potential. Within ACED, the exact KS correlation en- 
ergy can be written as 



E, = -J^""^l\ixldrj dr\Mr-r'y, 
{X\{r,r':iu) - xo{r , r' ; iu))} . 



(3) 



Here i'{r — r') = I/|r — r'| is the bare Coulomb interaction 
kernel, and xo is the KS independent-particle response 
function at imaginary frequencies iu, 



Xo(r,r-;..) = 2V"v^-(-)^"(-)^-(-')'^-(-') 



-f c.c. , 
(4) 
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where c.c. denotes "complex conjugate" and the prefactor 
2 acounts for the spin-degeneracy in closed-shell systems. 
In Eq. ([3]), xa is the density-density response function 
of the "intermediately" interacting many-electron system 
employing a scaled Coulomb potential I'x. We adhere to 
the commonly used notation of i,j ■■■ being occupied, 
i.e. hole KS states and a,b. .. being unoccupied or vir- 
tual, particle states. In principle, a Dyson-type integral 
equation^ has to be solved for xa, 

XA = XO + XO (I'A + /xc) XA , (5) 

with f^^ as the xc kernel, i.e. the functional derivative 
of the exchange-correlation potential with respect to the 
density. Within RPA, /xc = 0, i.e. using many-body 
terminology^ so-called vertex corrections are not in- 
cluded in the response function x or equivalcntly in the 
screening of the Coulomb interaction. Solving Eq. ([5]) for 
Xa with f^^ = corresponds to the diagrammatic resum- 
mation of ring grapha^i^ to infinite order. In passing 
we note that, working within RPA, Eq. ^ can be rear- 
ranged to 

XA = (1 -XoJ^a)"^ -Xo = [l + X0i^A + X0i^AX0'^A + - • •] -Xo, 

(6) 

reflecting the above mentioned summation of the 
(screened) Coulomb interaction up to infinite order 
in xoi^A- As will be seen later, Eq. (0) resem- 
bles the coupled-cluster amplitude equations where so- 
called particle-particle, particle-hole, and hole-hole lad- 
der terms have been removed (see Eq.ll7|). Starting from 
Eq. the A-integral is readily done and the final ex- 
pression for the RPA correlation energy reads 

/•OO J 

Ef"-^ = |^Tr{ln(l - xo(*")^) + Xo(m)^} • (7) 

B. From coupled-cluster theory to RPA and 
RPA+SOSEX 

From a DFT purist's point of view, the previously 
outlined ACFD terminology for the RPA is certainly 
the most consistent way to classify "RPA" as a corre- 
lation energy functional to the many-electron ground- 
state. An alternative formulation of the RPA may be 
motivated starting from many-body theory. Many-body 
or equivalcntly field-theoretical diagrammatic techniques 
originally developed in quantum electrodynamics and nu- 
clear physics^ have been applied to the homogeneous 
electron gas as well as finite systems like atoms and 
molecules for several decades already. For systems that 
are not strongly correlated, the most successful diagram- 
matic, partial summation technique (see Refs. |56| and 
[roh is the coupled cluster (CC) expansion of the many- 
electron wavefunction. The CC expansion to the ho- 
mogeneous electron gas has been applied by Freeman^ 
Kiimmcl, Liihrmann, and Zabolitzky;^ as well as Bishop 
and Liihrmann J ''^1^^ The same CC expansion techniques 



are indispensable ingredients for highly accurate molec- 
ular calculations. Here, pioneers have been Cizekf^iS 
Paldus et al.^ and Bartlett and Purvis^I to name a few. 
A more complete list of references may be found in the 
recent review article by Bartlett and Musiali^S 

The CC expansion relies on the ansatz for the many- 
electron wavefunction, j^*), 

m = e^m, (8) 

to generate the exact ground state from the ground state 
$) of the reference system commonly within the HF ap- 
proximation. Note that T may be represented by a sum 
of single, double, and higher-order excitation operators, 
generating in a similar way to configuration interaction 
(CI) techniques, singly, doubly substituted determinants 
based on the HE reference wavefunction |$). However, 
the CC expansion is distinct from CI by virtue of the 
exponential ansatz used in CC expansions (Eq. [5]) for the 
wavefunction j^*), with 

= 1 + f + if 2 + If 3 + . . . , (9) 

introducing so-called disconnected products of excita- 
tions responsible for the size-extensivity of the coupled 
cluster correlation energyjS 

In coupled-cluster doubles theory (CCD) the excita- 
tion operator corresponds to a double excitation operator 
only, where 

f = f2 with (10) 

The amplitudes f^j are obtained from solving a set of 
so-called doubles amplitude equations reading 

($^/|e-*i?e*|$) = 0. (12) 

Solving Eq. (|12p self-consistently for tf^ results in a re- 
summation of infinitely many diagrams of a certain type. 
Removing all terms from the above amplitude equation 
that do not correspond to so-called ring-diagrams defines 
the so-called ring-CCD. 

Recently, the equivalence between direct, 
i.e. "Coulomb term only" ring-CCD (drCCD) and 
RPA as considered by Freeman;^ reexamined by 
Griineis and Kresse^ and Scuseria et air- has been 
demonstrated. Scuseria et al. algebraically showed 
that the CCD approximation to the many-electron 
wavefunction contains the ring-approximation, i.e. the 
RPA to the ground-state correlation energy, but also 
includes selected higher-order exchange and ladder 
diagrams i^^iiSi^^ In other words, RPA equals drCCD 
and therefore corresponds to a subset of CCD diagrams. 

Within the framework of CC expansions, the RPA and 
RPA-f SOSEX correlation energies may be calculated us- 
ing drCCD amplitudes {t,°^} by employing the respective 



4 



equations 



6.36.37 



RPA 



E. 



E. 



RPA+SOSEX 



■ijab 



ah 



iajb ^ij 



(13) 
(14) 



ijab 



The matrices Biajb and Kiajb are of rank A^occ x -^virt, 
and they are defined by two-electron integrals Biajb = 
( ij I ab) and Kia.jb — {ij \ ab) — ( i j | 6a), respectively, 

{pq\rs) = //(/)*(x)(/),.(x)— ^— -(/)*(x')'?!)s(x')dxrfx', 



(15) 

with x={r,CT}. The amplitudes {ifj*} are obtained from 
solving a set of nonlinear Riccati equations, closely re- 
lated to to the time-dependent HF or more precisely the 
time-dependent Hartrec method)^ 

{ij\ab) + {Ec - ek)SacSikttj + {ic\ak)tfj 

+ t1^{e, ~ ek)5t,c5,k + tf^{ic\ak) (16) 
+ tt^{kl\cd)t'i^ = Q. 

The previous equation can be rewritten in a more com- 
pact form;^ 



B + AT + TA + TBT = 0, 



(17) 



with Aia,jb = (ca - ei)5ij6ab + {ib\aj), Bia.jb = (u|a6), 
and Tiajb = tfj, underlining the quadratic order in the 
amplitudes' matrix T. 

Freeman has evaluated the RPA correlation energy 
of the unpolarized electron gas for various electron 
densities^ using the drCCD equations and compared 
them to Hedin's RPA results (see Table II in Ref. [tI) 
following an approach suggested by Nozicres and Pines <^ 

Both agree to within the numerical accuracy employed 
in the calculations. Moreover, Freeman has gone beyond 
RPA via inclusion of the second-order screened exchange 
(SOSEX) diagram. He found that SOSEX reduces the 
correlation energy by about 30%. Monkhorst and Odd- 
ershede came to similar conclusions employing RPA and 
RPA+SOSEX to mctalhc hydrogen;^ and Griincis ob- 
served a similar reduction of the correlation energy for 
small atoms^ finding good agreement with highly accu- 
rate coupled cluster correlation energies only after inclu- 
sion of SOSEX. Finally we note that until recently the 
formulation of SOSEX within an AGED framework has 
not been entirely clear, but has lately been shown by 
Jansen et al.^ 



C. Single excitations and their renormalization 

As alluded to above, in most practical calcula- 
tions, RPA and SOSEX correlation energies are eval- 
uated using KS orbitals from local or semilocal den- 
sity functionalsj^ii^ or generalized KS orbitals^iii from 



range-separated density functionals. This way, both RPA 
and SOSEX can be interpreted in terms of many-body 
perturbation theory (MBPT) based on a (generalized) 
KS reference state, where only a selected type of di- 
agrams are summed up to infinite order. If one per- 
forms a simple Raylcigh-Schrodinger perturbation theory 
(RSPT) starting from an (approximate) KS-DFT refer- 
ence, and examines the perturbation series at second or- 
der, one can identify a term arising from single excita- 
tions (SE), that is not included in RPA or SOSEX corre- 
lation energies. In terms of single-particle orbitals, this 
term can be expressed as 



e: 



SE 



HF 



(18) 



cff • 

V IS 



where the self-consistent HF potential, and 

the effective single-particle potential that defines the non- 
interacting reference Hamiltonian giving rise to the 
single- particle orbitals \i) and \a) in the above expression. 
(See the supplemental material of Ref. [2^ for a detailed 
derivation.) As is obvious from Eq. E^^ trivially 

vanishes for the HF reference, i.e., when = . but 
is nonzero otherwise. It has been shown that adding 
this term to RPA improves the description of weak in- 
teractions significantly^^ Note that the choice of v'^^ in 
Eq. HH]) is slightly different in RSPT from that in the 
2nd-ordcr Gorling-Levy perturbation theory (GL2)i^ In 
the latter case, v''^ = t;Exx-OEP^ ^-^^-^ ^;Exx-oep 
the exact-exchange OEP— potential. The difference 
of the two perturbation theories lies in the choice of the 
adiabatic-connection path (A-integral) - in GL2 the elec- 
tron density is kept fixed along the path way and the 
perturbative Hamiltonian has a non-linear dependence 
on A, whereas in RSPT the A-dependence of the per- 
turbative Hamiltonian is linear, but the electron density 
varies along the A-integral. Eq. in RSPT is more 
efficient and practically useful in the sense that there is 
no need to solve the computationally intensive and some- 
times numerically problematic EXX-OEP equation and 
more flexible in the sense that it can be matched to any 
suitable reference state. The price one has to pay is that 
the theory, strictly speaking, is not KS-DFT formulated 
within the AGED framework. 

The SE contribution at second order as given by 
Eq. (fT5|) may become ill-behaved when the single-particle 
gap closes. To deal with this problem, in Ref. [2^ a se- 
quence of higher-order terms involving SE processes have 
been identified and summed up in the spirit of RPA. This 
leads to a "renormalized" SE (rSE) contribution to the 
correlation energy. 



^RSE ^ ^ . 



|(z|Az;|a) 



(i\/S.v\i) 



(a|Ai;|a)' 



(19) 



where Av 



,HF 



,,cff 



The additional term {i\/S.v\i) 



(a I Aw I a) in the denominator of Eq. ([T^ is negative def- 
inite, and prevents the possible divergence of the ex- 
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TABLE I: The list of methods used throughout this work and their acronyms, 
exact-exchange level is abbreviated by "EX". 



Note that the total energy at the 



(EX+RPA)@PBE 

HF+RPA@PBE 

(EX+RPA+SE)ePBE 

(EX+RPA+rSE)QPBE 

HF+(RPA+SOSEX)@PBE 

(EX+RPA+SOSEX)(aPBE 



EX and RPA evaluated with a PBE reference, i.e. PBE orbitals and eigenvalues 

HF total energy combined with RPA using a PBE reference 

EX and RPA augmented with SE using the PBE reference 

EX and RPA augmented with rSE using the PBE reference 

HF total energy combined with RPA+SOSEX using the PBE reference 

EX, RPA+SOSEX using the PBE reference 



(EX+RPA+SOSEX+rSE)@PBE EX, RPA+SOSEX, and rSE using the PBE reference 



pression even when the KS gap closes. The rSE cor- 
rection is therefore expected to have a more general 
applicability, while preserving the good performance of 
the 2nd-ordcr SE for wide-gap molecules and insulators. 
In deriving Eq. p^ . however, the "non-diagonal" ele- 
ments in the higher-order SE diagrams have been ne- 
glected for simplicity. Such an approach lacks invari- 
ance with respect to unitary transformations (orbital 
rotations) within the occupied and/or unoccupied sub- 
spaces. The orbital-rotation-invariance can be restored 
by including also the "non-diagonal" elements. This can 
be achieved by first scmi-diagonalizing the Fock Hamilto- 
nian / = h'^^ + v^^ — separately within the occupied 
and unoccupied subspaces of h'^^ and utilizing the resul- 
tant (so-called semi- canonical) orbitals and orbital ener- 
gies in Eq. (|18p . A detailed description of this procedure 
will be presented in a forthcoming paper. However, we 
emphasize that results presented in this work are based 
on Eq. ([T^ . but despite the lack of rotational invariance 
in the orbitals of this approach, numerical results arc only 
very little affected. 

As also demonstrated in Ref. [lH, the SE contributions 
to the correlation energy can be effectively accounted 
for to a large extent by replacing the non-sclf-consistent 
HF total energy computed using KS orbitals by its self- 
consistent counterpart. In this so-called hybrid- RPA 
scheme the RPA correlation energy is still evaluated us- 
ing KS orbitals, whereas the EX term is evaluated us- 
ing HF orbitals. The same strategy can be applied 
to "RPA+SOSEX" calculations. In this work, we will 
benchmark the influence of SE contributions on the per- 
formance of RPA and SO SEX both by explicitly includ- 
ing the (r)SE corrections and in terms of the hybrid 
scheme. 

As outlined in Ref. [2^ by Ren et al. , rendering the en- 
ergy functional stationary with respect to variations in 
the orbitals implies a zero correlation energy contribu- 
tion stemming from SEs. This is well known as Bril- 
louin's theorem. It will be demonstrated in this work, 
that SE effects represent a non-negligible contribution to 
the correlation energy and consequently affect results on 
thermochemistry and kinetics. In the field of quantum 
chemistry effects induced by SEs are known as orbital- 
relaxation effectsi^i^ Besides MBPT discussed above, 
the SE terms are present in the CC theory as well. In this 



context, Scuseria and Schaefer have shown that CCD em- 
ploying optimized-orbitals (see Ref. gives results very 
close to CCSD. On the other hand, optimizing orbitals 
for CCSD calculations docs not lead to significant im- 
provements in the wavefunction. In other words, changes 
in the correlation energy induced upon inclusion of SEs 
may be effectively incorporated by means of a unitary 
transformation, i.e. rotation of the orbitals, as given in 
Eq. (6) of Ref. Isi- 

We close this section by presenting Table HI which sum- 
marizes the acronyms of the various methods applied 
in this work. For the KS single-determinant reference 
wave function we use the Perdew, Burke, and Ernzerhof 
(PBE)21 generalized gradient approximation (GGA). We 
adopt the notation introduced by Ren et al. in Ref. ITtI . 
hence "@PBE" means "evaluated using PBE orbitals and 
orbital energies". This particular choice of orbitals is 
mainly driven by the following arguments: (i) PBE con- 
tains no empirically adjusted parameters, (ii) performs 
slightly better than LDA (see e.g. Ref. and (iii) 

it is computationally less expensive to calculate the or- 
bitals using semilocal functionals instead of e.g. hybrid 
functionals In addition, once one restricts the input or- 
bitals to KS orbitals, results have shown to be virtually 
identical to those obtained using PBE orbitalsi^i^ 



III. COMPUTATIONAL DETAILS 

Computational results of the present work are based 
on calculations using (i) the Vienna ah initio simula- 
tion package vasP;^"— (ii) a development version of the 
GAUSSiAr*22 suite of programs, and (iii) FHI-aims<^i2i 
All of the software packages used have the RPA and 
RPA+SOSEX functionals available since recentlyi^ii^iii 
VASP uses periodic boundary conditions and projector 
augmented plane waves as a basis set, which makes it ide- 
ally suited for extended, crystalline systems. GAUSSIAN is 
based on local, analytic Gaussian type (GT) basis func- 
tions using open boundary conditions and the linear com- 
bination of atomic orbitals to expand the molecular or- 
bitals. FHI-aims primarily uses numeric, atom-centered 
basis functions, but GT orbitals can be employed as well. 
In both cases, all the required integrals are evaluated nu- 
merically on an overlapping atom-centered gridi^ The 
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resolution-of-idcntity approximation is used to handle the 
four-centered Coulomb repulsion integrals and the KS re- 
sponse function (details of the implementation have been 
presented in Ref. lo^ . In this work GT orbitals are used 
in FHI-aims calculations to facilitate a direct comparison 
with GAUSSIAN and the extrapolations to the complete 
basis set (CBS) limit. 

In this work we present statistical errors for the 
G2-1 seti^^ii as well as for BH6,.^ HTBH38/04, and 
NHTBH38/04 sets of 38 hydrogen transfer and 38 non- 
hydrogen transfer barrier heights after Zhao et aZ4^i^ 
Results for the molecular test sets use a two-point ex- 
trapolation procedure on the correlation energies to at- 
tain the complete basis set (CBS) limit<^"— The chosen 
ansatz is motivated by an atomic partial wave expansion 
of the two-particle many-body wavefunctiou)^ 

^corr ~ ^corr ~^ ' i'^^) 

where the E^^^. are correlation energies corresponding to 
the cc-pVXZ basis sets. For G2-1, CBS calculations are 
based on Dunning's correlation-consistent cc-pVQZ and 
cc-pV5Z basis sets i^^i^°° Note that throughout this work 
CBS extrapolation will be denoted by, e.g., cc-pV(Q,5)Z. 

Moreover, G2-1 calculations employ the Boys-Bernardi 
counterpoise correctioniSi to correct for basis set su- 
perposition errors (BSSE) within a particular basis set. 
Therefore, we emphasize that the CBS procedure uses 
BSSE free correlation energies. In order to avoid inaccu- 
racies from numerical quadrature of xc energy contribu- 
tions, GAUSSIAN calculations use a grid of 400 radial shells 
and 770 angular points in each shell to converge the KS 
orbitals. GAUSSIAN employs a root-mean-square conver- 
gence criterion for the density matrix in the SCF iteration 
of 0.1 yuHartree, which implies an energy convergence no 
worse than at least 0.01 /^Hartree (GAUSSIAN keyword: 
SCF=tight). In FHI-aims the grid setting "tight" to- 
gether with "radial_multiplier=6" has been used to 
achieve convergence within one /xHartree. 

TABLE II: Matching radii Tc of the PAW potentials used in 
the present work. If the matching radii differ for specific quan- 
tum numbers, they are specified for each /-quantum number 
using subscripts. 





Valence 


Vc [a.u.] 




Valence 


Tc [a.u.] 


H 


Is 


1.0s l.lpd 


F 


2s2p 


1.1s 1.4pd 


Li 


ls2s 


1.2s I'Spd 


Mg 


2p3s 


2.0sd 1.6p 


B 


2s2p 


1.5s l.Tpd 


AI 


3s3p 


1.9spd 2.0/ 


C 


2s2p 


1.2s I'Spd 


Si 


3s3p 


1.5s 1.9pd 


N 


2s2p 


1.3s I'Spd 


p 


3s3p 


1.9sp 2.0d/ 





2s2p 


1.2s I'Spd 


CI 


3s3p 


1.7s l.9pdf 



Results on barrier heights in BH6, HTBH38/04, and 
NHTBH38/04 use a cc-pV(T,Q)Z CBS extrapolation of 
the correlation energies and do not employ counterpoise 
corrections. To test for the errors incurred, we com- 
pare with benchmark results obtained using RPA and 



TABLE III: Experimental lattice constants, a™f, extrapo- 
lated to K. Energy cutoffs for the one-electron wave func- 
tions -Epw as well as energy cutoffs for representing the over- 
lap charge densities employed in the calculation of the 
atomization energies of solids. The corresponding structures 
are denoted using the Strukturbericht symbols in parenthesis 
in the first column (A4=diamond, Bl=rock-salt, B3=zinc- 
blende). All energies and lattice constants in eV and A, re- 
spectively. 







Epw 


Ex 


C (A4) 


3.567=^ 


550 


400 


Si (A4) 


5.430=' 


450 


300 


SiC (B3) 


4.358"- 


550 


400 


BN (B3) 


3.607^ 


550 


400 


BP (B3) 


4.538'' 


450 


350 


AIN (B3) 


4.380^^ 


550 


400 


AlP (B3) 


5.460*^ 


450 


350 


LiH (Bl) 


4.064'' 


600 


450 


LiF (Bl) 


4.010'" 


600 


450 


LiCl (Bl) 


5.106"- 


600 


450 


MgO (Bl) 


4.207" 


600 


450 



"Ref. [To3, '"Ref. [lol, " Ref. [lol, Ref. Ho^ . 



RPA+SOSEX given in Ref. d. The statistical errors in 
barrier heights deviate from the aforementioned bench- 
mark values by at most 1 kJ/mol. Hence, the errors 
incurred using smaller basis sets are minute and conse- 
quently are not expected to bias the conclusions. 

The test set on atomization energies for crystalline 
solids includes 11 archetypal semiconductors and insu- 
lators. Specifically it comprises C, Si, SiC, BN, BP, 
AIN, AlP, LiH, LiF, LiCl, and MgO. The projector aug- 
mented wave (PAW) pseudopotentials (technical details 
in Tab. |ll| and kinetic energy cutoffs employed in the 
pres ent calculations are identical to the ones used in 
Ref. ll02l . Table IhII summarizes the lattice constants used 
in "post-RPA" calculations. Moreover, we specify plane 
wave cutoffs f or th e overlap charge densities described in 
Refs. m and Il02l . The SOSEX correlation energy was 
calculated using a (3 x 3 x 3) F-centered k mesh, except 
for BN and BP due to a slower /c-point convergence of 
the energy. For these systems a (4 x 4 x 4) mesh was 
used. RPA correlation energies are taken from the liter- 
ature (see Ref. Issh . In VASP, atoms are calculated using 
a supercell approach. The dimension of the supercells 
has been chosen as (9 x 9 x 9) in size. To reduce the 
computational cost of the "RPA+SOSEX" calculations 
for isolated atoms, natural orbitals obtained using sec- 
ond order pertu rbati on theory have been employed. As 
outlined in Ref. Il07l natural orbitals substantially im- 
prove convergence of the correlation energy with respect 
to the number of virtual orbitals. 

To assess the codes used in this work, we com- 
pare numerical results obtained using the "RPA" and 
"RPA+SOSEX" implementations of GAUSSIAN and FHI- 
aims. Table IIVI shows correlation energies for the He 
atom obtained using the cc-pV5Z basis set. In order 
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TABLE IV: Benchmark calculations for atomic He using 
FHI-aims and GAUSSIAN and a cc-pV5Z GT orbital basis 
set. Results are given in Hartree atomic units. 



He / cc-pV5Z 


GAUSSIAN 


FHI-aims 


HP 


-2.86162468 


-2.86162483 


MP2 


-0.03640606 


-0.03640651 


RPA@HF 


-0.06524488 


-0.06524570 


(RPA+SOSEX)@HF 


-0.03262244 


-0.03262285 



to avoid errors caused by numerical integration, we de- 
cided to use (restricted, i.e. spin-unpolarizcd) HF or- 
bitals and eigenvalues for the calculation of RPA and 
RPA-I-SOSEX. The agreement found is close to per- 
fect. Differences between results are within a sub-micro- 
Hartree error margin. In passing we mention that FHI- 
aims employs the resolution-of-identity (RI) techniquep^ 
which (i) reduces the computational workload signifi- 
cantly and (ii), as shown in Tab. IIVI docs not sacrifice 
accuracy. For the molecular test sets, we always cross- 
check the "RPA" and "RPA-hSOSEX" results obtained 
with the GAUSSIAN suite of program and FHI-aims to 
make sure that the results presented in this work are not 
affected by the actual implementations. "SE" and "rSE" 
have so far only been implemented in FHI-aims and we 
use these results throughout. 



IV. RESULTS AND DISCUSSION 

Central findings of this work arc summarized in Tab.lVl 
presenting binding energies in molecules (G2-1) and 
solids, HT activation energies or barrier heights (BH6, 
HTBH38) as well as NHT barrier heights (NHTBH38). 
Whenever results are compared to experiment or to the 
best theoretical estimates, we use mean error (ME) and 
mean unsigned error (MUE) as statistical measures to as- 
sess the accuracy of individual methods employed. Note 
that the experimental reference values are corrected for 
zero poin t effects and are taken from the liter atur e (G2-1: 
Ref. llOSi : atomization energies in solids: Ref . Il09|; barrier 
heights in BH6, HT/NHTBH38: Refs. il Elandlii). 
Reaction energies, as presented in Table I VII are calcu- 
lated from barrier heights in HTBH38 and NHTBH38, 
respectively. 



A. Atomization energies of small molecules and 
solids 

The notorious underbinding of (EX+RPA)@PBE in 
molecules and solids has already been demonstrated in 
many studies .^^'^'i^ii^'ii'^ii^ Table |V] presents MEs and 
MUEs in binding (atomization) energies obtained us- 
ing (EX-f-RPA)@PBE for insulating sohds (see Sec. Hill 
as well as for the molecules contained in the G2-1 set. 



8 - (EX+RPA)@PBE 




FIG. 1: Mean unsigned relative errors (MURE) in the atom- 
ization energies of 55 small molecules contained in G2-1 (full 
bars) and 11 insulating solids (squared bars) obtained using 
four of the RPA-based methods presented in this work. 



(EX+RPA)@PBE 
HF+RPA@PBE 
{EX+RPA+SE)@PBE 
(EX+RPA+rSE)(3)PBE 




FIG. 2: Mean unsigned relative errors (MURE, given in %) 
in G2-1 for all RPA-based methods presented in this work. 
Atomization energies use counterpoise correction and corre- 
lation energies are CBS extrapolated using cc-pV(Q,5)Z. 

On average, (EX-|-RPA)@PBE underbinds solids com- 
pared to experiment by — 67kJ/mol (see also Ref.HH) and 
molecules by — 43kJ/mol. We repeat that experimental 
binding energies arc corrected for zero-point eff ects and 
are taken from the literature (for G2-1, see Ref. [M for 
the test set on solids see Ref. Il09l and references therein) . 

Following the suggestion of Ren et al.,— effects in- 
curred by replacing the EX@PBE reference energy by the 
HF total energy have been checked for both molecules 
and sohds. Indeed, HF-|-RPA@PBE improves binding 
energies of molecules and solids by almost 50% compared 
to (EX4-RPA)@PBE. Fig. [T] presents mean unsigned rel- 
ative errors (MURE) in molecular (full bars) as well as 
solid state (squared bars) binding energies. Overall dif- 
ferences in MUREs arc rather small. For the commonly 
applied (EX-hRPA)@PBE method, the MURE is found 
to be approximately 6%. Using HF at the exact ex- 
change level reduces the MURE by more than 1%. It 
appears that the aforementioned improvements are less 
pronounced at the relative scale, and the error reduction 
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TABLE V: Mean errors (ME) and mean unsigned errors (MUE) in atomization or binding energies of 11 solids (see Tab. 
and 55 molecules (G2-1), in the barrier heights comprised in BH6, in HTBH38/04 (hydrogen transfer barriers), as well as in 
NHTBH38/04 (non- hydrogen transfer barriers). Results are given in kj/mol. 



Method 


Sohds 


G2-1 


BH6 


HTBH38 


NHTBH38 




ME MUE 


ME 


MUE 


ME MUE 


ME MUE 


ME MUE 


(EX+RPA)@PBE 


-67.5 


67.5 


-42.7 


42.8=" 


1.2 7.5=' 


-0.8 


7.1 


-10.5 12.1 


HF+RPA@PBE 


-34.7 


36.7 


-25.3 


30.3 


-25.5 25.5 


-36.8 


36.8 


-48.5 48.5 


(EX+RPA+SE)@PBE 






-14.2 


22.9 


-23.8 23.8 


-52.7 


52.7 


-50.6 51.9 


(EX+RPA+rSE) @PBE 






-26.2 


27.4 


-14.8 16.3 


-18.0 


21.7 


-31.4 31.4 


(EX+RPA+SOSEX)@PBE 


-27.0 


27.0 


-20.3 


23.0" 


17.6 17.6^ 


22.2 


22.2 


13.4 15.5 


HF+(RPA+SOSEX)@PBE 


5.8 


17.4 


-2.9 


13.0 


-9.2 9.2 


-13.8 


13.8 


-24.7 25.5 


(EX+RPA+SOSEX+rSE)@PBE 






-4.0 


13.9 


3.1 3.7 


3.6 


5.4 


-6.3 17.6 



See Ref. [^. Note t hat differences in the MUE of G2-1 are due to the different values for the experimental dissociation 
energies (see Ref. llOSi ). 



is apparently bigger for solids than for molecules. 

The explicit inclusion of the SE contribution to the cor- 
relation energy "SE@PBE" obtained using Eq. (|T8)) has 
been evaluated for molecules only. Adding "SEQPBE" 
to (EX+RPA)@PBE results in an ME of approximately 
— 14kJ/mol (see Tab. W} and an MUE of approximately 
23kJ/mol, clearly outperforming HF+RPA@PBE. Rel- 
ative unsigned errors in G2-1 collected in Fig. [2] further 
corroborate the improvements of (EX-|-RPA+SE)@PBE 
over HF-|-RPA(S)PBE. Overall, these results confirm the 
findings presented by Ren et al. in Ref. [2^ However, 
"Renormalization" of the SE contributions, as required 
for systems with a small one-electron band gap in PBE 
(see activation energies discussed in Sec. I IV Bp brings the 
atomization energies in the G2-1 test set back into almost 
perfect agreement with HF+RPA@PBE. Therefore, the 
good agreement with experiment for the G2-1 test set on 
the level of (EX-hRPA-^SE)@PBE is most likely to some 
extent fortuitous. 

As extensively discussed in Refs. 113 and H, 
the (RPA-f-SOSEX) correlation energy, here denoted 
as "(RPA+SOSEX)@PBE," represents a correction 
to (EX-|-RPA)@PBE rectifying the one-electron self- 
interaction error contained in "RPA@PBE" due to exclu- 
sion principle violating diagrams <^ Results for G2-1 ob- 
tained using (EX+RPA+SOSEX)@PBE are taken from 
Ref. i and included in Tab. El The ME in G2-1 obtained 
using (EX-HRPA+SOSEX)@PBE is approximately equal 
to — 20kJ/mol. For solids, the ME error reduces to 
-27kJ/mol. Compared to (EX+RPA)@PBE, this rep- 
resents substantial improvements of approximately 50% 
for atomization energies. 

Given that both SOSEX and rSE, or alternatively 
replacing EXQPBE by HF, alleviate the tendency of 
(EX-|-RPA)@PBE to underbind, both schemes are ex- 
pected to work cooperatively for the atomization ener- 
gies of small molecules. Indeed, replacing "EXQPBE" 
in (EX-hRPA-HSOSEX)@PBE by the HF total en- 
ergy yields excellent results, with a slight underbind- 
ing for molecules (ME = — 2.9kJ/mol), and a slight 
overbinding for solids (ME=:5.8 kJ/mol). Again, the 



HF-1-(RPA+S0SEX)@PBE (ME = -2.9kJ/mol) and 
the (EX+RPA+SOSEX+rSE)@PBE methods (ME = 
— 4.0kJ/mol) perform almost on par for molecules. 

In summary, single excitation diagrams improve 
(EX-fRPA)@PBE atomization energies of small 
molecules at virtually zero additional computational 
cost. However, as we will see below, this method fails 
when the one-electron band gaps in PBE become small. 
The better founded rSE does not perform equally well for 
atomization energies when comined with RPA. Combined 
with RPA+SOSEX it yields impressive atomization en- 
ergies that arc also in almost perfect agreement with 
the "hybrid variants" e.g. the (self-consistent) HF total 
energy together with "(RPA+SOSEX)@PBE" . Overall 
this indicates that (EX-|-RPA+SOSEX+rSE)@PBE or 
HF+(RPA+SOSEX)@PBE are the methods of choice 
for atomization energies. 



B. Activation energies in HTBH38 and NHTBH38 
chemical reactions 



Transition State 



Reactants 



Activation i = Forward 
Energy "^Barrier (Vp 



Reaction ' 
Energy \ 
(AE) • 



Products 



Reverse 
Barrier (Vf ) 



Reaction Coordinate 
FIG. 3: Schematic of activation and reaction energies. 
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The ability to accurately describe the topology of mul- 
tidimensional potential energy surfaces spanned by the 
internal molecular degrees of freedom, i.e. the reaction 
coordinates, in the course of a chemical reaction, is 
central to first principles electronic structure methods. 
Calculating the energy difference between reactants and 
transition states (see Fig. [3]) is a stringent test for the ac- 
curacy of density functionals. As mentioned in Sec. IIIIl 
the HTBH38 and NHTBH38 test sets established by 
Truhlar and coworkcrs^SiiS will be used here to test the 
RPA-bascd functionals considered. 

Our findings on barrier heights, i.e. activation energies 
(Fig. [3]), are summarized in Table IVl MEs and MUEs are 
calculated with respect to the best theoretical estimates 
currently available for HT and NHT barrier heights given 
in Rcfs.l48landl49l respectively. Furthermore the MUREs 
in HT barriers [panel (a)] and NHT barriers [panel (b)] 
are depicted in Fig. H) Note that legends given in Fig. |4] 
follow the color code used in Fig. [21 To establish a con- 
nection to Ref. i. Tab. |V] also shows MEs and MUEs 
for the BH6 test set,— which has been introduced as 
a computationally less intensive, but statistically rep- 
resentative subset of HT/NHTBH38. However, we do 
not present a detailed discussion on BH6 here, but stress 
that errors in BH6 essentially follow the trends found for 
HT/NHTBH38. 




FIG. 4: Panel a) shows mean unsigned relative errors 
(MURE) in HT barrier heights of HTBH38. Panel b) shows 
MUREs in NHTBH38 for the RPA-based methods presented 
in this work. Energies use a cc-pV(T,Q)Z extrapolation and 
the frozen core approximation in calculated correlation ener- 
gies. 

One of the main findings of this work is the 
astonishingly good performance of the conventional 
(EX+RPA)@PBE scheme for activation energies. To be 
more specific, (EX-f RPA)@PBE performs significantly 
better for the transfer of hydrogen atoms than for reac- 



tions involving heavier atoms. For HTBH38, the ME ob- 
tained using (EX+RPA)@PBE amounts to -0.8kJ/mol 
and the associated MUE amounts to 7.1kJ/mol. These 
error margins are similar to those of some of the range- 
separated, generalized KS-DFT functionals like e.g. LC- 
cjPBE.iiiS The latter performs very well for chemical 
reaction barriers (see also Sec. lIVPp . However, for 
(EX-i-RPA)@PBE, the MUE increases by more than 
50% when elements heavier than H, like e.g. F or CI, 
are transferred. The MUE in NHTBH38 obtained us- 
ing (EX+RPA)@PBE amounts to 12.1kJ/mol together 
with a rather distinct underestimation of the barriers by 
-10.5kJ/mol (ME). 

On a relative scale, the MURE for HT reactions ob- 
tained using (EX+RPA)@PBE (see Fig. g]) amounts to 
approximately 20%, but increases to a value approxi- 
mately twice as large for NHT reactions [panel (b)] . Note 
that reaction 7 in NHTBH38 has a barrier height of only 
— 1.42kJ/mol. For this reaction MUREs are extraordi- 
narily large leading to an increase, which is seven to eight 
times as large as the corresponding value in HT reactions. 
The statistics would be drastically biased by such a case 
being very likely compensated by significantly extending 
the test set. Therefore, we decided to exclude reaction 
number 7 from the test set when calculating the MURE 
in NHTBH38. 

Both HF-|-RPA@PBE and (EX+RPA+SE)@PBE 
show a strong underestimation of barriers with maxi- 
mal errors as large as 50kJ/mol. As mentioned above, 
the reason for this behavior has been attributed to small 
HOMO-LUMO differences found for some of the transi- 
tion state structures, which arc not correctly described by 
the simple (EX+RPA-hSE)@PBE scheme. Indeed, the 
renormalization of SE alleviates the problem, and the 
corresponding ME and MUE in HTBH38 obtained us- 
ing (EX+RPA+rSE)@PBE are improved by almost 60% 
compared to (EX+RPA+SE)@PBE. Note that numer- 
ical results given in Tab. |V] nicely reflect the trend in- 
duced by incorporation of SE effects in the correlation 
energy contribution, i.e. it partially takes care of the 
lack of self-consistency in the EX@PBE energy. However, 
the rSE corrects for the strong underestimation of barri- 
ers seen in HF+RPA@PBE and (EX+RPA-f SE)@PBE, 
but qualitatively reflects the same trend compared to 
(EX+RPA)@PBE. 

The performance of (EX-t-RPA-fSOSEX)@PBE for 
barrier heights has already been tested by Paier et 
al. for the BH6 test seti^ This work extends the 
findings of Ref. |^ by discriminating HT and NHT 
reactions. (EX-HRPA+SOSEX)@PBE is less accu- 
rate for HT barriers than (EX-HRPA)@PBE as indi- 
cated by an MUE of about 22kJ/mol compared to 
7kJ/mol. Quantitatively, (EX+RPA+SOSEX)@PBE 
on average overestimates barrier heights for HTBH38 
by the aforementioned 22kJ/mol. This is in per- 
fect agreement with the errors found for the BH6 test 
set^ On the other hand, (EX+RPA+SOSEX)@PBE 
performs substantially better for NHT barrier heights. 
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where ME and MUE are found to be close to the 
ones obtained usmg (EX+RPA)@PBE. On average, 
(EX+RPA+SOSEX)@PBE overestimates NHT barriers 
by approximately 13kJ/mol, whereas (EX+RPA)@PBE 
underestimates them by llkJ/mol. As shown in 
Fig. m the MURE in NHT barriers obtained using 
(EX+RPA+SOSEX)@PBE amounts to 34% [panel (b) 
in Fig. U shghtly outperforming (EX+RPA)@PBE by 
approximately 3%. 

Incorporation of SE effects into 
(EX+RPA+SOSEX)@PBE in the hybrid fashion, 
i.e. HF+(RPA+SOSEX)@PBE, leads to very different 
results when applied to HT and NHT reactions, respec- 
tively. HF+(RPA+SOSEX)@PBE improves HT reaction 
barrier heights, whereas NHT barrier heights deteriorate 
appreciably compared to (EX+RPA+SOSEX)@PBE, 
ending up with an overall underestimation of barrier 
heights. 

The situation becomes noticeably better, for 
both HT and NHT barrier heights, upon combina- 
tion of explicitly computed renormalized SE with 
(EX-|-RPA+SOSEX)@PBE. Barrier heights obtained 
using (EX-l-RPA-fSOSEX+rSE)@PBE are of similar 
quality as "conventional" (EX-f RPA)@PBE, although 
the unsigned error in NHT test set is slightly larger. 
(EX+RPA+SOSEX+rSE)@PBE overestimates HT 
barriers by approximately 3.6kJ/niol, but reduces the 
ME in NHT barriers (ME = -6.3kJ/mol) compared to 
(EX-hRPA)@PBE. 

To summarize this section, SOSEX and rSE tend to 
overestimate and underestimate reaction barrier heights, 
respectively. Thus it appears advantageous to com- 
bine the correction schemes in order to achieve a par- 
tial error compensation. This mechanism works most 
efficiently for HT reactions and somewhat less so for 
NHT reactions. Taking the excellent performance 
of (EX+RPA+SOSEX+rSE)@PBE for binding energies 
(see previous Section) into account, this functional of- 
fers the most balanced description in terms of binding 
energies as well as activation energies. 



C. Reaction energies in HTBH38 and NHTBH38 

As shown in Fig. [31 knowing both forward (V^) and 
reverse (V^) reaction barrier heights, corresponding re- 
action energies AE are readily calculated using 

AE = vf-Vf. (21) 

Note that 17 out of the 38 reactions contained in 
HTBH38 lead to a nonzero AE, whereas NHTBH38 com- 
prises 13 reactions with a forward barrier different from 
the reverse barrier. The corresponding MEs and MUEs 
of the RPA-based functionals are compiled in Tab. IVIl 
and the MUREs are depicted in Fig. [S] 

Similar to the trends found for atomization energies, 
HT reaction energies are significantly improved upon in- 
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FIG. 5: Panel a) shows mean unsigned relative errors 
(MURE) in HT reaction energies of HTBH38. Panel b) shows 
MURE for the reaction energies of NHTBH38. Energies use 
a cc-pV(T,Q)Z extrapolation and the frozen core approxima- 
tion to the correlation energies. Color code of legends follows 
Fig.gl 



elusion of (SOSEX)@PBE to (EX+RPA)@PBE as re- 
fiected in the MUEs. For (EX+RPA)@PBE the MUE in 
HT reactions amounts to 18.2kJ/mol and drops down 
to 4.6kJ/mol employing (EX+RPA+SOSEX)@PBE. 
Hence, it appears that eliminating the one-electron self- 
correlation error contained in RPAQPBE is beneficial 
for HT reaction energies. This is not entirely surpris- 
ing, since the aforementioned error will be largest for 
breaking and creating covalent hydrogen bonds. For re- 
actions involving heavier atoms, as exemplified by the re- 
action energies in NHTBH38, the correction due to (SO- 
SEX) @PBE appears to perform less favorably. This can 
be seen by inspection of Fig. [5] presenting MUREs in HT 
[panel (a)] as well as NHT reaction energies [panel (b)]. 
For (EX-hRPA)@PBE the MUE in NHTBH38 amounts 
to 9.7kJ/mol, which is rather low, whereas for NHT reac- 
tion energies obtained using (EX+RPA+SOSEX)@PBE, 
the MUE increases to 20.5kJ/mol. 

Concerning effects due to SE@PBE and rSEOPBE 
to (EX-f-RPA)@PBE, no significant improvement of 
HT reaction energies over (EX-|-RPA)@PBE has been 
found. The MEs and MUEs given in Tab. |Vl] for 
(EX-hRPA-hSE)@PBE (ME = -3 kJ/mol; MUE = 
16.9 kJ/mol) and (EX+RPA-hrSE)@PBE (ME = -2.9 
kJ/mol; MUE ~ 17 kJ/mol) are essentially unaltered 



11 



TABLE VI: Mean errors and mean unsigned errors [kj/mol] in the reaction energies 
obtained using calculated barrier heights of the HTBH38/04 hydrogen transfer as well 
as NHTBH38/04 non- hydrogen transfer barrier heights. 



HTBH38 NHTBH38 



Method 


ME 


MUE 


ME 


MUE 


(EX+RPA)@PBE 


-3.2 


18.2 


-7.8 


9.7 


HF+RPA®PBE 


2.2 


12.3 


-1.6 


14.4 


(EX+RPA+SE)@PBE 


-3.0 


16.9 


9.4 


24.6 


(EX+RPA+rSE) @PBE 


-2.9 


17.0 


-1.2 


11.8 


(EX+RPA+SOSEX)@PBE 


2.7 


4.6 


-18.4 


20.5 


HF+(RPA+SOSEX)@PBE 


2.8 


4.1 


-12.2 


15.7 


(EX+RPA+SOSEX+rSE)@PBE 


3.0 


4.4 


-11.9 


15.5 



compared to (EX+RPA)@PBE. In contrast to HT, the 
rSE correction helps to improve the NHTBH38 reac- 
tion energies and alleviates the overestimation found for 
simple (EX+RPA+SE)@PBE drastically (ME = -1.2 
kJ/mol compared to 9.4 kJ/mol). The associated MUE 
as well as MURE decrease by approximately 50%. 

We now turn to a discussion of results obtained 
using the "hybrid variants" , which employ the HE 
energy as the reference energy on the EX level. 
Specifically for (EX-|-RPA)@PBE, HT reaction en- 
ergies are substantially improved upon replacement 
of EXQPBE through HE. As can be seen from 
Tab. IVIl the MUE is reduced by approximately 
6kJ/mol, which translate into an improvement of the 
MURE by approximately 50%. HT reaction energies 
obtained using (EX+RPA+SOSEX)@PBE, which 
are fairly accurate, are hardly affected by changing 
to the corresponding hybrid scheme. Employing 
HF-H(RPA+SOSEX)@PBE, however, the MUE in NHT 
reaction energies is reduced by 5kJ/mol. In addition, 
the ME amounts to — 12kJ/mol, which compares very 
favorably to the ME of — 18kJ/mol obtained using 
(EX-|-RPA+SOSEX)@PBE. In terms of performance, 
the combined scheme (EX+RPA+SOSEX+rSE)@PBE 
is on par with HF-f(RPA+SOSEX)@PBE for 
both HTBH38 and NHTBH38 reaction energies. 
(EX-FRPA+SOSEX+rSE)@PBE has two apparent 
favorable features: (i) it substantially improves HT 
reaction energies obtained using (EX+RPA)@PBE, 
and (ii) it performs approximately similarly well for 
all of the test sets investigated in this work. In 
other words, the overall variation in error margins 
for atomization energies, barrier heights, and reaction 
energies is smallest for (EX+RPA+SOSEX+rSE)@PBE 
lending the functional robustness. Among the func- 
tionals discussed in this work, (EX+RPA)@PBE 
performs best for NHT reaction energies. Neverthe- 
less, (EX+RPA+SOSEX+rSE)@PBE performs only 
slightly worse, but given the better HT reaction barrier 
heights and the significantly better reaction energies in 
HTBH38, (EX+RPA+SOSEX-HrSE)@PBE is among 
the RPA-based functionals tested in this work, the 
functional of broadest applicability. 



D. Comparing RPA to semilocal and hybrid 
functionals 

To close the discussion on the performance of the 
RPA- and RPA-|-SOSEX-based functionals, we briefly 
compare molecular atomization and activation energies 
to results obtained using commonly applied semilo- 
cal as well as HF/DFT hybrid functionals. To ren- 
der direct comparisons easier. Table IVIII repeats MUEs 
for G2-1, BH6, HTBH38, and NHTBH38 for three 
of the RPA-based functionals, which perform best, 
namely (EX+RPA)@PBE, HF-h(RPA-|-SOSEX)@PBE, 
and (EX+RPA+SOSEX-hrSE)@PBE. The above men- 
tioned statistical errors are compared to PBE-GGA, 
BLYP-GGA"3,ii4 as well as the PBEOiii^ and 
B3LYP-ii^ global hybrid functionals. In addition, we also 
present results obtained using the above mentioned LC- 
cjPBE range-separated hybrid functional^^ LC-wPBE 
mixes a fraction of EX at the long-range of the Coulomb 
interaction (for definitions, see Ref. IllOl ). but uses only 
one parameter (0.40 bohr~^) for defining a universal 
range separation. It is remarkable that LC-wPBE de- 
scribes reaction barriers and atomization energies ex- 
tremely accurately representing a landmark among hy- 
brids for thermochemistry and kinetics. Admittedly, for 
extended systems admixture of EX on the long-range is 
detrimental and leads to e.g. strongly overestimated band 
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gapSliii 

Returning to RPA, activation energies obtained using 
(EX-fRPA)@PBE are de facto on par with LC-wPBE 
(Tab. Trends for GGA and global hybrid func- 

tionals like PBEO or B3LYP are rather general, hence 
other GGA-type and global hybrid functionals pe rform 
quite similarly (for other functionals, see e.g. Ref. Ill2h . 
Although, HE+(RPA-|-SOSEX)@PBE does not perform 
as well as (EX+RPA)@PBE for activation energies of 
non-hydrogen transfer reactions (corresponding MUE 
is almost twice as large), it performs certainly better 
than PBE and BLYP. HF-|-(RPA+SOSEX)@PBE is only 
slightly outperformed by B3LYP for the aforementioned 
activation energies in NHTBH38. According to this syn- 
opsis, (EX-FRPA+SOSEX+rSE)@PBE certainly shows 
the most balanced description of molecular binding and 



12 



TABLE VII: Comparing the three best-performing functionals presented in Tab. |V] 
to widely used semilocal and HF/DFT hybrid functionals. Mean unsigned errors in 
individual test sets are given in kj/mol. 



Method 


G2-1 


BH6 


HTBII38 


NIITBH38 


(EX+RPA)@PBE 


42. 8=^ 


7.5" 


7.1 


12.1 


HF+(RPA+SOSEX)@PBE 


13.0 


9.2 


13.8 


25.5 


(EX+RPA+SOSEX+rSE)@PBE 


13.9 


3.7 


5.4 


17.6 


PBE 


36.0*' 


38.9" 


39.0"* 


33.9"* 


BLYP 


19.7*= 


32.6" 


31.5"* 


36.4"* 


PBEO 


14.6*= 


19.2" 


17.7"* 


14.1"* 


B3LYP 


10.0*= 


19.7" 


17.7"* 


18.2"* 


LC-cjPBE 


15.6° 




5.4" 


10.0" 



See Ref. Note that differences in the MUE of G2-1 a re du e to the different 
valu es fo r the experimental dissociation energies (see Ref. Il08l l . 

^ Ref. fll3 

" Ref. m 

<* Ref. 4^^ 

" Ref. I lid . Note that the MUE given here for G2 refers to G2-2 comprising 148 
molecules. The MUE for G2-1 will be lower. 



barrier heights. It performs similarly well as hybrid 
functionals in terms of atomization energies, outperforms 
both GGA and hybrid functionals in terms of hydrogen- 
transfer barrier heights, and performs equivalently well 
for non-hydrogen barrier heights as aforementioned hy- 
brids do. 

Although, this work is not devoted to weak, van-dcr- 
Waals-type of interactions, it should be emphasized that 
all of the RPA-based functionals presented here do in- 
clude them seamlessly as already mentioned in the intro- 
duction. It is well known that neither GGA nor hybrid 
functionals do show the correct van der Waals asymp- 
tote. 



V. CONCLUSIONS 

In summary, we have reported an extensive assess- 
ment of several exact-exchange based post-KS density 
functionals involving RPA correlation energies and be- 
yond. Correlation energies have been assessed for solids 
as well as for small molecules. Specifically we cal- 
culated atomization energies of solids and molecules 
using (EX-hRPA)@PBE, (EX-|-RPA+SOSEX)@PBE as 
well as HF+RPA@PBE and HF+(RPA+SOSEX)@PBE, 
where the latter approach gives binding energies im- 
proved by approximately 50% compared to "conven- 
tional" (EX-|-RPA)@PBE. Furthermore, we investigated 
the performance of individual functionals for chemi- 
cal reaction barrier heights or activation energies em- 



ploying large and well established test sets. Gener- 
ally, we found that it is rather difficult to system- 
atically improve on (EX+RPA)@PBE reaction bar- 
rier heights, although modest improvements using 
(EX-hRPA-fSOSEX-hrSE)@PBE were achieved for HT 
barriers. Importantly, the favorable impact of the corre- 
lation energy contribution stemming from SE effects on 
binding energies does not translate to reaction barriers. 
This has been explained by divergent correlation energy 
contributions in systems with small HOMO-LUMO gaps. 
Therefore, application of "SE" to systems with small one- 
electron band gaps is not possible, but a renormaliza- 
tion of "SE" helps to alleviate the problem. Surprisingly, 
(EX-f RPA)@PBE yields reaction energies of high accu- 
racy for reactions involving non-hydrogen atoms. Good 
and robust performance of a novel RPA-based functional 
(EX-hRPA-fSOSEX-hrSE)@PBE is a central finding of 
this work. It improves on binding or atomization en- 
ergies compared to (EX-|-RPA)@PBE, improves on HT 
barrier heights as well as reaction energies. 
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